rm(list = ls())
require(stats4)
require(foreign)
require(extRemes)



data21=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december6noon_type7clean.csv", header=TRUE)
data22=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december15_type7clean.csv", header=TRUE)
data23=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_january26_type7(clean).csv", header=TRUE)

dataframe21=data.frame(data21)
dataframe22=data.frame(data22)
dataframe23=data.frame(data23)


#some descriptive stats and graphs
exp2frame=rbind(dataframe21,dataframe22,dataframe23)
exp2frame=subset(exp2frame,uid>800)
exp2frame$correct <-(exp2frame$state==4&(exp2frame$response==2|exp2frame$response==7|exp2frame$response==8|exp2frame$response==9))|(exp2frame$state==3&(exp2frame$response==1|exp2frame$response==6|exp2frame$response==3|exp2frame$response==5))
exp2frame$incentive <- (exp2frame$qid==5)*5+(exp2frame$qid==1)*40+(exp2frame$qid==6)*70+(exp2frame$qid==7)*95
accuracyexp2=c(sum((exp2frame$incentive==5)*exp2frame$correct)/sum(exp2frame$incentive==5),sum((exp2frame$incentive==40)*exp2frame$correct)/sum(exp2frame$incentive==40),sum((exp2frame$incentive==70)*exp2frame$correct)/sum(exp2frame$incentive==70),sum((exp2frame$incentive==95)*exp2frame$correct)/sum(exp2frame$incentive==95))
incentivec=c(5,40,70,95)
uidlistexp2=unique(exp2frame$uid)

#check N
#**
length(uidlistexp2)
#**

exp2totcor=c(sum(exp2frame$correct*(exp2frame$incentive==5)),sum(exp2frame$correct*(exp2frame$incentive==40)),sum(exp2frame$correct*(exp2frame$incentive==70)),sum(exp2frame$correct*(exp2frame$incentive==95)))
exp2totinc=c(sum((1-exp2frame$correct)*(exp2frame$incentive==5)),sum((1-exp2frame$correct)*(exp2frame$incentive==40)),sum((1-exp2frame$correct)*(exp2frame$incentive==70)),sum((1-exp2frame$correct)*(exp2frame$incentive==95)))


#find the multinomial constant
multinom=function(vec){
if(length(vec)==2){
return(sum(log(1:max(1,sum(vec))))-sum(log(c(1:max(1,vec[1]),1:max(1,vec[2])))))
}
if(length(vec)==3){
return(sum(log(1:max(1,sum(vec))))-sum(log(c(1:max(1,vec[1]),1:max(1,vec[2]),1:max(1,vec[3])))))
}
}



exp2multinomvec=vector("numeric")
for(i in 1:4){
exp2multinomvec[i]=multinom(c(exp2totcor[i],exp2totinc[i]))
}
exp2multinomconst=sum((exp2multinomvec))

######################
#linear mutual information

exp2acculin=function(kappaall, kappaloc){
accuvec=vector('numeric')
for(i in 1:4){
cost=function(accu){
return((kappaall+kappaloc)*(accu*log(accu)+(1-accu)*log(1-accu)))
}
negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-cost(accu)))
}
optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

exp2likelihoodlin=function(kappaall,kappaloc){
accuvec=exp2acculin(kappaall,kappaloc)
like=exp2multinomconst+sum(exp2totcor*log(accuvec)+exp2totinc*log(1-accuvec))
return(-like)
}

######################
#linear mutual information no neighborhood

exp2acculinnoneigh=function(kappaall){
accuvec=vector('numeric')
for(i in 1:4){
cost=function(accu){
return((kappaall)*(accu*log(accu)+(1-accu)*log(1-accu)))
}
negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-cost(accu)))
}
optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

exp2likelihoodlinnoneigh=function(kappaall){
accuvec=exp2acculinnoneigh(kappaall)
like=exp2multinomconst+sum(exp2totcor*log(accuvec)+exp2totinc*log(1-accuvec))
return(-like)
}


################################
#General Entropy Model

exp2accugen=function(kappaall,kappaloc,rho){
accuvec=vector('numeric')
for(i in 1:4){

cost=function(accu){
if(rho==1){
totcost=(kappaall+kappaloc)*(accu*log(accu)+(1-accu)*log(1-accu))
}
if(rho==2){
totcost=(kappaall+kappaloc)*(-1/2*(log(accu)+log(1-accu))-2*log(2))
}
if(rho!=1&rho!=2){
totcost=(kappaall+kappaloc)*(2^(1-rho)/((rho-1)*(rho-2))*(accu*accu^(1-rho)+(1-accu)*(1-accu)^(1-rho))-1/((rho-1)*(rho-2))-log(2))
}
return(totcost)
}


negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-cost(accu)))
}

optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

exp2likelihoodgen=function(kappaall,kappaloc,rho){
accuvec=exp2accugen(kappaall,kappaloc,rho)
like=exp2multinomconst+sum(exp2totcor*log(accuvec)+exp2totinc*log(1-accuvec))
return(-like)
}

####################################
#Power mutual information

exp2accupow=function(kappaall,kappaloc,rho){
accuvec=vector('numeric')
for(i in 1:4){

cost=function(accu){
totcost=(accu*log(accu)+(1-accu)*log(1-accu)-2*0.5*log(0.5))
if(rho>0){
return(totcost^rho)
}else{
return(-1*totcost^rho)
}
}

negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-(kappaall+kappaloc)*cost(accu)))
}

optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

exp2likelihoodpow=function(kappaall,kappaloc,rho){
accuvec=exp2accupow(kappaall,kappaloc,rho)
like=exp2multinomconst+sum(exp2totcor*log(accuvec)+exp2totinc*log(1-accuvec))
return(-like)
}

######################################################################################################################################################################################

#import data
symmetrydata=read.dta("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Expt1.dta")
symmetryframe=data.frame(id=symmetrydata$var1,user_id=symmetrydata$var2,qn=symmetrydata$var3,qid=symmetrydata$var4,order=symmetrydata$var5,bid=symmetrydata$var6,chosen_act=symmetrydata$var9, true_state=symmetrydata$var10)

#reward value
r=10 


#data frame for letters
symframeletters=subset(symmetryframe,user_id<1700)
symframeletters$correct=(symframeletters$true_state<=13 & symframeletters$chosen_act==1)+(symframeletters$true_state>13 & symframeletters$chosen_act==2)
unique(symframeletters$qid)

#data frame for balls
symframeballs=subset(symmetryframe,user_id>1700)
symframeballs$true_state=symframeballs$true_state+6*(symframeballs$qid==3)+4*(symframeballs$qid==4)+2*(symframeballs$qid==5)
symframeballs$correct=(symframeballs$true_state<=10 & symframeballs$chosen_act==1)+(symframeballs$true_state>10 & symframeballs$chosen_act==2)
unique(symframeballs$qid)

#Variables to be stored
acculin=matrix(rep(0,80),nrow=4,ncol=20)
accupow=matrix(rep(0,80),nrow=4,ncol=20)
accugen=matrix(rep(0,80),nrow=4,ncol=20)
accunoneigh=matrix(rep(0,80),nrow=4,ncol=20)
dataaccu=matrix(rep(0,80),nrow=4,ncol=20)

#Qids for each treatment type
lettersqid=c(3,5,7,8)
ballsqid=c(3,4,5,6)

#number of states used by different treatments
Nvec=c(8,12,16,20)

#States for different treatments
states=list()
states[[1]]=c(7:14)
states[[2]]=c(5:16)
states[[3]]=c(3:18)
states[[4]]=c(1:20)

correct=list()
incorrect=list()

#{*} Make correct vecs (needs changing between treatment types)
for(i in 1:4){
correct[[i]]=vector('numeric')
incorrect[[i]]=vector('numeric')
for(j in 1:Nvec[i]){
correct[[i]][j]=sum(symframeballs$correct*(symframeballs$qid==ballsqid[i])*(symframeballs$true_state==states[[i]][j]))
incorrect[[i]][j]=sum((1-symframeballs$correct)*(symframeballs$qid==ballsqid[i])*(symframeballs$true_state==states[[i]][j]))
}
}



exp4multinomlist=list()
exp4multinomvec=vector("numeric")
for(i in 1:4){
exp4multinomlist[[i]]=vector("numeric")
for(j in 1:Nvec[i]){
exp4multinomlist[[i]][j]=multinom(c(correct[[i]][j],incorrect[[i]][j]))
}
exp4multinomvec[i]=sum((exp4multinomlist[[i]]))
}

exp4multinomconst=sum(exp4multinomvec)

######################
#Linear Mutual Information Costs

exp4accuracylin=function(kappaall,kappaloc,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
overallcost=kappaall*sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))
neighborcost=vector('numeric')
for(j in 1:(N/2-1)){
neighborcost[j]=overallcost=kappaall*sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))
}
neighborcost[N/2]=kappaloc*(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))

}
neighborcost[N/2]=kappaloc*(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))
return(overallcost+sum(neighborcost))
}

#Optimization target
negpayoff=function(accu){
return(-1*(r*sum(accu/(N/2))-cost(accu)))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2),negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

exp4likelihoodlin=function(kappaall,kappaloc){
accuvec1=exp4accuracylin(kappaall,kappaloc,1)
accuvec2=exp4accuracylin(kappaall,kappaloc,2)
accuvec3=exp4accuracylin(kappaall,kappaloc,3)
accuvec4=exp4accuracylin(kappaall,kappaloc,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=exp4multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}
 
#################################################

#power costs

exp4accuracypow=function(kappaall,kappaloc,rho,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
overallcost=sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))-N*1/N*log(1/N)
neighborcost=vector('numeric')
for(j in 1:(N/2-1)){
neighborcost[j]=1/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))-2*0.5*log(0.5)
}
neighborcost[N/2]=(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))-2*0.5*log(0.5)
return(kappaall*overallcost^rho+kappaloc*sum(neighborcost^rho))
}

#Optimization target
negpayoff=function(accu){
return(-1*(r*sum(accu/(N/2))-cost(accu)))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2)-c(1:(N/2))/200,negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

exp4likelihoodpow=function(kappaall,kappaloc,rho){
accuvec1=exp4accuracypow(kappaall,kappaloc,rho,1)
accuvec2=exp4accuracypow(kappaall,kappaloc,rho,2)
accuvec3=exp4accuracypow(kappaall,kappaloc,rho,3)
accuvec4=exp4accuracypow(kappaall,kappaloc,rho,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=exp4multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}

#################################################
#Gen entropy approach
exp4accuracygen=function(kappaall,kappaloc,q,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
neighborcost=vector('numeric')

if(q==1){

overallcost=kappaall*sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))
}
neighborcost[N/2]=kappaloc*(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))

}else{

if(q==2){
overallcost=kappaall*(-1/N*sum(log(accu/(N/2))+log((1-accu)/(N/2)))-2*log(N))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc*(-1/2*((accu[j]+accu[j+1])/2*(log(accu[j]/(accu[j]+accu[j+1]))+log(accu[j+1]/(accu[j]+accu[j+1])))+(2-accu[j]-accu[j+1])/2*(log((1-accu[j])/(2-accu[j]-accu[j+1]))+log((1-accu[j+1])/(2-accu[j]-accu[j+1]) ) )) -2*log(2))
}
neighborcost[N/2]=kappaloc*(-1/2*(log(accu[N/2])+log(1-accu[N/2]))-2*log(2))


}else{
overallcost=kappaall*(N^(1-q)/((q-2)*(q-1))*(sum((accu/(N/2))^(2-q)+((1-accu)/(N/2))^(2-q)))-1/((q-2)*(q-1))-log(N))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc*(2^(1-q)/((q-2)*(q-1))*((accu[j]+accu[j+1])/2*((accu[j]/(accu[j]+accu[j+1]))^(2-q)+(accu[j+1]/(accu[j]+accu[j+1]))^(2-q))+(2-accu[j]-accu[j+1])/2*(((1-accu[j])/(2-accu[j]-accu[j+1]))^(2-q)+((1-accu[j+1])/(2-accu[j]-accu[j+1]))^(2-q)))-1/((q-2)*(q-1))-log(2))

}
neighborcost[N/2]=kappaloc*(2^(1-q)/((q-2)*(q-1))*((accu[N/2])^(2-q)+(1-accu[N/2])^(2-q))-1/((q-2)*(q-1))-log(2))

}

}

return(overallcost+sum(neighborcost))
}

#Optimization target
negpayoff=function(accu){
return(-1*(r*sum(accu/(N/2))-cost(accu)))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2),negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

exp4likelihoodgen=function(kappaall,kappaloc,q){
accuvec1=exp4accuracygen(kappaall,kappaloc,q,1)
accuvec2=exp4accuracygen(kappaall,kappaloc,q,2)
accuvec3=exp4accuracygen(kappaall,kappaloc,q,3)
accuvec4=exp4accuracygen(kappaall,kappaloc,q,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=exp4multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}

######################
#Linear Mutual Information Costs no neighborhoods

exp4accuracynoneigh=function(kappaall,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Information cost function
cost=function(accu){
overallcost=kappaall*(accu*log(accu/(N/2))+(1-accu)*log((1-accu)/(N/2)))
return(overallcost)
}

#Optimization target
negpayoff=function(accuscalar){
return(-1*(r*sum(rep(accuscalar,N/2)/(N/2))-cost(accuscalar)))
}

#Find optimal accuracies
optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec=rep(optimout$minimum,N/2)
return(accuvec)
}

exp4likelihoodnoneigh=function(kappaall){
accuvec1=exp4accuracynoneigh(kappaall,1)
accuvec2=exp4accuracynoneigh(kappaall,2)
accuvec3=exp4accuracynoneigh(kappaall,3)
accuvec4=exp4accuracynoneigh(kappaall,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=exp4multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}
 
######################################################################################################################################################################################


#Multi-experiment numbers**********************************
likelihoodlin=function(kappaall,kappaloc){
return(exp2likelihoodlin(kappaall,kappaloc)+exp4likelihoodlin(kappaall,kappaloc))
}

modellin=mle(likelihoodlin,start=list(kappaall=3,kappaloc=24),method="L-BFGS-B",lower=c(0.001,0),upper=c(100,100))
modellin

likelihoodpow=function(kappaall,kappaloc,rho){
return(exp2likelihoodpow(kappaall,kappaloc,rho)+exp4likelihoodpow(kappaall,kappaloc,rho))
}

modelpow=mle(likelihoodpow,start=list(kappaall=14.1,kappaloc=88.9,rho=2.01),method="L-BFGS-B",lower=c(0.001,0,0.001),upper=c(1000,1000,1000))
modelpow

likelihoodgen=function(kappaall,kappaloc,rho){
return(exp2likelihoodgen(kappaall,kappaloc,rho)+exp4likelihoodgen(kappaall,kappaloc,rho))
}

modelgen=mle(likelihoodgen,start=list(kappaall=0.02,kappaloc=1.8,rho=8.4),method="L-BFGS-B",lower=c(0.001,0,0.001),upper=c(100,100,100))
modelgen

likelihoodnoneigh=function(kappaall){
return(exp2likelihoodlinnoneigh(kappaall)+exp4likelihoodnoneigh(kappaall))
}
modelnoneigh=mle(likelihoodnoneigh,start=list(kappaall=4),method="L-BFGS-B",lower=c(0.001),upper=c(100))
modelnoneigh


#Some LR tests
noneightest=lr.test(-logLik(modelnoneigh),-logLik(modellin),alpha=0.05,df=1)
gentest=lr.test(-logLik(modellin),-logLik(modelgen),alpha=0.05,df=1)
powtest=lr.test(-logLik(modellin),-logLik(modelpow),alpha=0.05,df=1)


########################################################################################

#General Entropy Model

exp2accugennoneigh=function(kappaall,rho){
accuvec=vector('numeric')
for(i in 1:4){

cost=function(accu){
if(rho==1){
totcost=(kappaall)*(accu*log(accu)+(1-accu)*log(1-accu))
}
if(rho==2){
totcost=(kappaall)*(-1/2*(log(accu)+log(1-accu))-2*log(2))
}
if(rho!=1&rho!=2){
totcost=(kappaall)*(2^(1-rho)/((rho-1)*(rho-2))*(accu*accu^(1-rho)+(1-accu)*(1-accu)^(1-rho))-1/((rho-1)*(rho-2))-log(2))
}
return(totcost)
}


negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-cost(accu)))
}

optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

exp2likelihoodgennoneigh=function(kappaall,rho){
accuvec=exp2accugennoneigh(kappaall,rho)
like=exp2multinomconst+sum(exp2totcor*log(accuvec)+exp2totinc*log(1-accuvec))
return(-like)
}

####################################
#Power mutual information

exp2accupownoneigh=function(kappaall,rho){
accuvec=vector('numeric')
for(i in 1:4){

cost=function(accu){
totcost=(accu*log(accu)+(1-accu)*log(1-accu)-2*0.5*log(0.5))
if(rho>0){
return(totcost^rho)
}else{
return(-1*totcost^rho)
}
}

negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-(kappaall)*cost(accu)))
}

optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

exp2likelihoodpownoneigh=function(kappaall,rho){
accuvec=exp2accupownoneigh(kappaall,rho)
like=exp2multinomconst+sum(exp2totcor*log(accuvec)+exp2totinc*log(1-accuvec))
return(-like)
}

############Single Experiment Numbers************************

modellinexp2only=mle(exp2likelihoodlinnoneigh,start=list(kappaall=24),method="L-BFGS-B",lower=c(0.001),upper=c(100))
summary(modellinexp2only)

modelpowexp2only=mle(exp2likelihoodpownoneigh,start=list(kappaall=7728,rho=4),method="L-BFGS-B",lower=c(0.001,0.001),upper=c(100000,1000))
summary(modelpowexp2only)

modelgenexp2only=mle(exp2likelihoodgennoneigh,start=list(kappaall=2,rho=6),method="L-BFGS-B",lower=c(0.001,0,0.001),upper=c(100,100,100))
summary(modelgenexp2only)

modellinexp4only=mle(exp4likelihoodlin,start=list(kappaall=5,kappaloc=5),method="L-BFGS-B",lower=c(0.001,0),upper=c(100,100))
summary(modellinexp4only)

modelpowexp4only=mle(exp4likelihoodpow,start=list(kappaall=5.1,kappaloc=4.69,rho=0.93),method="L-BFGS-B",lower=c(0.001,0,0.001),upper=c(1000,1000,1000))
summary(modelpowexp4only)

modelgenexp4only=mle(exp4likelihoodgen,start=list(kappaall=5.3,kappaloc=5,q=1.05),method="L-BFGS-B",lower=c(0.001,0,0.001),upper=c(100,100,100))
summary(modelgenexp4only)

modelnoneighexp4only=mle(exp4likelihoodnoneigh,start=list(kappaall=5),method="L-BFGS-B",lower=c(0.001),upper=c(100))
summary(modelnoneighexp4only)



#Likelihood Ratio Tests
1-pchisq(832, 1)

1-pchisq(39, 2)

1-pchisq(72,2)

1-pchisq(46,2)
